{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.8 Discontinuous Galerkin Methods\n",
    "\n",
    "* Use discontinuous finite element spaces to solve PDEs. \n",
    "* Allows upwind-stabilization for convection-dominated problems\n",
    "* Requires additional jump terms for consistency \n",
    "\n",
    "Interior penalty DG form for $-\\Delta u$:\n",
    "\n",
    "$$\n",
    "\\DeclareMathOperator{\\Div}{div}\n",
    "A(u,v) = \\sum_T \\int_T \\nabla u \\nabla v\n",
    "-  \\sum_F \\int_F \\{ n \\nabla u \\} [v] \n",
    "-  \\sum_F \\int_F \\{ n \\nabla v \\} [u] \n",
    "+ \\frac{\\alpha p^2}{h} \\sum_F \\int_F [u][v]\n",
    "$$\n",
    "\n",
    "with jump-term over facets:\n",
    "$$\n",
    "[u] = u_{left} - u_{right}\n",
    "$$\n",
    "\n",
    "and averaging operator\n",
    "$$\n",
    "\\{ n \\nabla u \\} = \\tfrac{1}{2} (n_{left} \\nabla u_{left} + n_{left} \\nabla u_{right})\n",
    "$$\n",
    "\n",
    "DG form for $\\Div (b u)$, where $b$ is the given wind:\n",
    "\n",
    "$$\n",
    "B(u,v) = -\\sum_T b u \\nabla v + \\sum_F \\int_F b\\cdot n   u^{upwind} v \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import netgen.gui\n",
    "from netgen.geom2d import unit_square\n",
    "from ngsolve import *\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The space is responsible for allocating the matrix graph. Tell it that it should reserve entries for the coupling terms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "order=4\n",
    "fes = L2(mesh, order=order, dgjumps=True)\n",
    "u,v = fes.TnT()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Every facet has a master element. The value from the other element is referred to via the\n",
    "`Other()` operator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "jump_u = u-u.Other()\n",
    "jump_v = v-v.Other()\n",
    "n = specialcf.normal(2)\n",
    "mean_dudn = 0.5*n * (grad(u)+grad(u.Other()))\n",
    "mean_dvdn = 0.5*n * (grad(v)+grad(v.Other()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Integrals on facets are computed by setting `skeleton=True`. \n",
    "* `dx(skeleton=True)` iterates over all internal faces\n",
    "* `ds(skeleton=True)` iterates over all boundary faces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 4\n",
    "h = specialcf.mesh_size\n",
    "a = BilinearForm(fes)\n",
    "diffusion = grad(u)*grad(v) * dx \\\n",
    "    +alpha*order**2/h*jump_u*jump_v * dx(skeleton=True) \\\n",
    "    +(-mean_dudn*jump_v-mean_dvdn*jump_u) * dx(skeleton=True) \\\n",
    "    +alpha*order**2/h*u*v * ds(skeleton=True) \\\n",
    "    +(-n*grad(u)*v-n*grad(v)*u)* ds(skeleton=True)\n",
    "\n",
    "a += diffusion\n",
    "a.Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = LinearForm(fes)\n",
    "f += SymbolicLFI(1*v)\n",
    "f.Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes, name=\"uDG\")\n",
    "gfu.vec.data = a.mat.Inverse() * f.vec\n",
    "Draw (gfu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DG requires a lot of additional matrix entries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a nze: 18900\n",
      "a2 nze: 5400\n"
     ]
    }
   ],
   "source": [
    "print (\"a nze:\", a.mat.nze)\n",
    "fes2 = L2(mesh, order=order)\n",
    "a2 = BilinearForm(fes2)\n",
    "a2 += SymbolicBFI(u*v)\n",
    "a2.Assemble()\n",
    "print (\"a2 nze:\", a2.mat.nze)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we are solving a convection-diffusion problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 4\n",
    "h = specialcf.mesh_size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `IfPos` checks whether the first argument is positive. Then it returns the second one, else the third one. This is used to define the upwind flux. The check is performed in every integration-point on the skeleton:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = CoefficientFunction( (20,5) )\n",
    "uup = IfPos(b*n, u, u.Other())\n",
    "\n",
    "convection = -b * u * grad(v)*dx + b*n*uup*jump_v * dx(skeleton=True)\n",
    "\n",
    "acd = BilinearForm(fes)\n",
    "acd += diffusion + convection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "acd.Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes)\n",
    "gfu.vec.data = acd.mat.Inverse(freedofs=fes.FreeDofs(),inverse=\"umfpack\") * f.vec\n",
    "Draw (gfu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hybrid Discontinuous Galerkin methods\n",
    "use additionally the *hybrid* facet variable on the skeleton:\n",
    "\n",
    "$$\n",
    "\\DeclareMathOperator{\\Div}{div}\n",
    "A(u,\\widehat u; v, \\widehat v) = \n",
    "  \\sum_T \\int_T \\nabla u \\nabla v\n",
    "- \\sum_T \\int_{\\partial T} n \\nabla u (v-\\widehat v)\n",
    "- \\sum_T \\int_{\\partial T} n \\nabla v (u-\\widehat u)\n",
    "+ \\frac{\\alpha p^2}{h} \\sum_F \\int_F (u-\\widehat u)(v-\\widehat v)\n",
    "$$\n",
    "\n",
    "the jump-term is now replaced by the difference $u - \\widehat u$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No additional matrix entries across elements are produced. Dirichlet boundary conditions are set as usual to the facet variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "order=4\n",
    "V = L2(mesh, order=order)\n",
    "F = FacetFESpace(mesh, order=order, dirichlet=\"bottom|left|right|top\")\n",
    "fes = FESpace([V,F])\n",
    "u,uhat = fes.TrialFunction()\n",
    "v,vhat = fes.TestFunction()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the jump is the difference between element-term and facet-term:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "jump_u = u-uhat\n",
    "jump_v = v-vhat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A non-zero elements: 4650\n"
     ]
    }
   ],
   "source": [
    "alpha = 4\n",
    "condense = True\n",
    "h = specialcf.mesh_size\n",
    "n = specialcf.normal(mesh.dim)\n",
    "\n",
    "a = BilinearForm(fes, condense=condense)\n",
    "dS = dx(element_boundary=True)\n",
    "a += grad(u)*grad(v)*dx + \\\n",
    "    alpha*order**2/h*jump_u*jump_v*dS + \\\n",
    "    (-grad(u)*n*jump_v - grad(v)*n*jump_u)*dS\n",
    "\n",
    "b = CoefficientFunction( (20,1) )\n",
    "uup = IfPos(b*n, u, uhat)\n",
    "a += -b * u * grad(v)*dx + b*n*uup*jump_v *dS\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(fes)\n",
    "f += SymbolicLFI(1*v)\n",
    "f.Assemble()\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "\n",
    "print (\"A non-zero elements:\", a.mat.nze)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not condense:\n",
    "    inv = a.mat.Inverse(fes.FreeDofs(), \"umfpack\")\n",
    "    gfu.vec.data = inv * f.vec\n",
    "else:\n",
    "    f.vec.data += a.harmonic_extension_trans * f.vec \n",
    "    \n",
    "    inv = a.mat.Inverse(fes.FreeDofs(True), \"umfpack\")\n",
    "    gfu.vec.data = inv * f.vec\n",
    "    \n",
    "    gfu.vec.data += a.harmonic_extension * gfu.vec\n",
    "    gfu.vec.data += a.inner_solve * f.vec\n",
    "\n",
    "Draw (gfu.components[0], mesh, \"u-HDG\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Projected Jumps\n",
    "The polynomial order of the facet-space can be reduced by one by inserting an $L_2$-projection into the penalty term. This can be implemented by keeping the highest order basis functions of the facet space discontinuous, see [Lehrenfeld-Schöberl: High order exactly divergence-free Hybrid Discontinuous Galerkin Methods for unsteady incompressible flows](https://www.sciencedirect.com/science/article/pii/S004578251630264X?via%3Dihub):\n",
    "\n",
    "$$\n",
    "\\DeclareMathOperator{\\Div}{div}\n",
    "A(u,\\widehat u; v, \\widehat v) = \n",
    "  \\sum_T \\int_T \\nabla u \\nabla v\n",
    "- \\sum_T \\int_{\\partial T} n \\nabla u (v-\\widehat v)\n",
    "- \\sum_T \\int_{\\partial T} n \\nabla v (u-\\widehat u)\n",
    "+ \\frac{\\alpha p^2}{h} \\sum_F \\int_F (\\Pi u-\\widehat u)(\\Pi v-\\widehat v)\n",
    "$$\n",
    "\n",
    "This is achieved by setting the flag `highest_order_dc=True` in the `FacetFESpace`. \n",
    "The number of matrix entries is reduced, while the order of convergence is preserved.\n",
    "This trick does not work in the case of significant convection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A non-zero elements: 2976\n"
     ]
    }
   ],
   "source": [
    "order=4\n",
    "V = L2(mesh, order=order)\n",
    "F = FacetFESpace(mesh, order=order, dirichlet=\"bottom|left|right|top\", \\\n",
    "                          highest_order_dc=True)\n",
    "fes = FESpace([V,F])\n",
    "u,uhat = fes.TrialFunction()\n",
    "v,vhat = fes.TestFunction()\n",
    "\n",
    "jump_u = u-uhat\n",
    "jump_v = v-vhat\n",
    "\n",
    "alpha = 2\n",
    "h = specialcf.mesh_size\n",
    "n = specialcf.normal(mesh.dim)\n",
    "\n",
    "a = BilinearForm(fes, condense=True)\n",
    "dS = dx(element_boundary=True)\n",
    "a += grad(u)*grad(v)*dx + \\\n",
    "    alpha*(order+1)**2/h*jump_u*jump_v*dS + \\\n",
    "    (-grad(u)*n*jump_v - grad(v)*n*jump_u)*dS\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(fes)\n",
    "f += SymbolicLFI(1*v)\n",
    "f.Assemble()\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "\n",
    "f.vec.data += a.harmonic_extension_trans * f.vec \n",
    "    \n",
    "inv = a.mat.Inverse(fes.FreeDofs(True), \"sparsecholesky\")\n",
    "gfu.vec.data = inv * f.vec\n",
    "    \n",
    "gfu.vec.data += a.harmonic_extension * gfu.vec\n",
    "gfu.vec.data += a.inner_solve * f.vec\n",
    "\n",
    "Draw (gfu.components[0], mesh, \"u-HDG\")\n",
    "\n",
    "print (\"A non-zero elements:\", a.mat.nze)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remarks on sparsity pattern in NGSolve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Remark 1: The sparsity pattern is set up a-priorily\n",
    "* The sparsity pattern of a sparse matrix in NGSolve is independent of its entries (it's set up a-priorily). \n",
    "* We can have \"nonzero\" entries that have the value 0\n",
    "\n",
    "Below we show the reserved memory for the sparse matrix and the (numerically) non-zero entries in this sparse matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes2 = L2(mesh, order=order, dgjumps=True)\n",
    "u,v=fes2.TnT()\n",
    "a3 = BilinearForm(fes2)\n",
    "a3 += u*v*dx + (u+u.Other())*v*dx(skeleton=True)\n",
    "a3.Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 1200x1200 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.sparse as sp\n",
    "import matplotlib.pylab as plt\n",
    "plt.rcParams['figure.figsize'] = (12, 12)\n",
    "A = sp.csr_matrix(a3.mat.CSR())\n",
    "fig = plt.figure(); ax1 = fig.add_subplot(121); ax2 = fig.add_subplot(122)\n",
    "ax1.set_xlabel(\"numerically non-zero\"); ax1.spy(A)\n",
    "ax2.set_xlabel(\"reserved entries (potentially non-zero)\"); ax2.spy(A,precision=-1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Remark 2: Sparsity pattern with and without `dgjumps=True` is different"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAADFCAYAAABXVUnDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAGZhJREFUeJzt3X/0XHV95/HnSyCBbRCEpPAt6iYu2XXRaoQvP9xajxVRfrgNnLI09qxGZTfuLmy1u7aGylp6gD2xVVk9rXTDgiHWFhDrkhraGpGzHPfIj4TGkCDQaOJqNhCikJKlggnv/eN+vjB8mfl+586dO3N/vB7nfM935t47n/nMzPt+5nPvfc/no4jAzMya72XjroCZmY2GG3wzs5Zwg29m1hJu8M3MWsINvplZS7jBNzNri4gY6x9wNvAwsB1YWaCcncADwGZgY1p2DLAB+Lv0/xWzlHEDsAfY2rGsaxmAgM+lem8BTs5R5hXArlTXzcC5HesuS2U+DLyrR5mvAu4EHgS2AR8uWtcZyixa18OBe4HvpHJ/Py1fBNyTHn8zMCctn5vub0/rF+Yocw2wo6OuS/J8Vo5tx3a/da1rXJce9LME4SHA94DXAHPSCz2pwE4xf9qyPyDtaMBK4JOzlPFW4ORpAdy1DOBc4K/Sm34GcE+OMq8APtpl25PSezA3Bc73gEO6bDcx9eECRwKPpMcOXNcZyixaVwHz0u3DUrCfAdwCLEvL/wT49+n2fwD+JN1eBtyco8w1wIVdtu/rs3JsO7b7rWtd43rcp3ROA7ZHxPcj4lngJmDpEMtfCtyYbt8InD/TxhFxF/CTPstYCqyNzN3A0ZIm+ixzpvreFBHPRMQOsm/u07qUuTsi7k+3nwK+C5xQpK4zlFm0rhER+9Pdw9JfAG8Hbu1R16nXcCtwpiT1WeZMdZ31sxoyx/ZLn6sxsV3XuB53g38C8MOO+z9i5g9iJgF8XdImSSvSsuMiYne6/Shw3ADl9iqjaN0vlbRF0g2SXjFomZIWAm8i6w0Mpa7TyixcV0mHSNpMdvi/gazH9GREHOjy2OfLTev3AcfOVmZETNX16lTXayTNHeT1D4lju+GxXce4HneDP0xviYiTgXOASyS9tXNlZMdAhcaRGEYZybXAPwGWALuBTw9SiKR5wFeAj0TE33euG7SuXcosXNeIOBgRS4BXkvWUXpu3jNnKlPR6svOurwVOJTvn+7Giz1MRju0OVYntOsb1uBv8XWQXVKa8Mi3LLSJ2pf97gK+SfQCPTR3ipP97Bii6VxkD1z0iHksf7HPAdbxwuNh3mZIOIwveL0XEXwyjrt3KHEZdO173k2QXz95Mdvh5aJfHPl9uWn8U8OM+yjw7HbpHRDwDfKFIXYfAsd2S2K5TXI+7wb8PWCxpkaQ5ZBcz1uUtRNLPSTpy6jbwTmBrKmt52mw5cNsAdexVxjrgfcqcAezrOOScrb6d59kuSHWdKnOZpLmSFgGLya7aT3+8gOuB70bEZ4ZR115lDqGuCyQdnW4fAZxFdg71TuDCHnWdeg0XAt9MPbrZynyoo0EQ2bnTzroO9FkV4NhucGzXNq6j5GyF2f7IrjQ/Qnb+6+MDlvEasqvqU+lMH0/LjwXuIEvl+gZwzCzl/DnZod3PyM6HXdyrDLIr43+c6v0AMJmjzC+mx2xJH9pEx/YfT2U+DJzTo8y3kB3SbqEjpaxIXWcos2hd3wD8bXr8VuATHZ/ZvWQXxL4MzE3LD0/3t6f1r8lR5jdTXbcCf8oLGQ99fVaObcd2v3Wta1wrPdDMzBpu3Kd0zMxsRNzgm5m1hBt8M7OWcINvZtYSpTX4ks6W9LCk7ZJWlvU8ZqPkuLY6K6XBl3QIWbrQOWQDEb1H0kkzbL+i17qC9Rh6uS6znWWmcnPFdVl1aXOZZZXbljLL6uHnHTiqlMAoqVyX2c4yYbAB0ery+upSZlnltqLMshr8cQxWZVY2x7XVWik/vJJ0IdkYEP8m3X8vcHpEXNqxzQrSt9XLjnj5KYce9fNDr4fl94snHNVz3eOPP86CBQuG+nyjKHPTpk17I6Lwk/QT12m5Y9uGamq/LBrbh86+yUBmHdQnIlYDqwHmTiyOieX/raSqWB57Z1inaevnz5vDxsvPKrlGxUn6wZCK6muwKse2DdvUfte5D86fNwc2vTNXbJfV4D8/cBTZDrEM+I2SnsvGZO/+Z1m4cn3hcuryxYHj2ipk7/5ncz+mlAY/Ig5IuhT4G7Kp3m6IiG1lPJfVX94vjp2rziuxNr05rq1q5hx/4il5ti+rh09E3A7cXlb51l55vhw0+ya5OK6tzvxLW2u0IH8vyKyp3OCbmbVEaad0qmCmc72TV20Y6KKHmVldNbrBn8mwskIWrVw/lJmfzczK1toGf1h25MwY8ReEmY2LG/wRy/sF0csw8t/NrF0a3eD3ahTHlcc9THleg69XmBk0vMG3TN7rFf6CaIepToNPM7aHG3x7CV/QbpdhnWbspg6dh52rzmtNrLrBt9L4grbl7TyM+trU/HlzgPyxWtdraG7wrTL8BWHDMIprdL2eo+pfBG7wrbb6/YLQJ9+9qeSqmAGjTwjJG9utbPC7fQuLcs9lmlVNk7PYrLtCDb6kncBTwEHgQERMSjoGuBlYCOwELoqIJ4pVs3w+NWCdmhTbZlOGMXjar0TEkoiYTPdXAndExGLgjnTfrI4c29YoZYyWuRS4Md2+ETi/hOcwGwfHttVa0XP4AXxdUgD/Pc3leVxE7E7rHwWO6/bAzome5xx/Ys8nmH4+sepXwa0xSo9ts1Er2uC/JSJ2Sfp5YIOkhzpXRkSkHeYlpk/0XLAeZsPm2LbGKdTgR8Su9H+PpK8CpwGPSZqIiN2SJoA9RZ5j0cr1I8uecfaOTRlFbFvzFD0DUXaG1MDn8CX9nKQjp24D7wS2AuuA5Wmz5cBtRSo47u7RuJ/fRm9UsW02akV6+McBX5U0Vc6fRcRfS7oPuEXSxcAPgIuKV9NspBzb1kgDN/gR8X3gjV2W/xg4s0ilzMbJsW1NVbtf2tZ1DAszm51/5VuuMvLwzcysgmrXwx8HjzlibZL3aNn7QX24h29m1hJu8M3MWqIxDb7GXQEzs4przDn8Xr+GdfaOmVmmMT18MzObWS16+AtXrq9kJoCzd8yqqepj2oyLe/hmZi3hBt/MrCUa3+A7e8fMLDPrOXxJNwDvBvZExOvTsq6TOSsbXvCzwLnA08D7I+L+cqreH2fvWC9VjO1u544dq/VR9XP//Vy0XQP8EbC2Y9nUZM6rJK1M9z8GnAMsTn+nA9em/63SbQedP28OGy8/awy1sRmsoWKxvXDlesdKhSxauX7gOTGqOHnSrKd0IuIu4CfTFveazHkpsDYydwNHp5mBWm/v/mfHXQWbpqqx7VipjiITIFVx8qRBz+H3msz5BOCHHdv9KC17CUkrJG2UtHHAOpiVwbFtjVX4om1EBAN8mUXE6oiYjIjJonUwK4Nj25pm0Ab/sanD2WmTOe8CXtWx3SvTsspx9o71UPvYNutl0F/aTk3mvIoXT+a8DrhU0k1kF7T2dRweV4qzd6yHysV23swPx7D10k9a5p8DbwPmS/oR8HtkO0O3yZxvJ0tb206WuvaBEupcW912xKqncTWZY3s4usV1FTNU8mjql+asDX5EvKfHqpdM5pzOeV5StFJmo+DYLk8VM1SsBb+0NTOzjBt8M7OWcIM/jbN3zKypajEe/ig5e8fMmqo2Df7kVRsaOb6IJ1GxphpGJ8njCg1XbU7peHwRs/bxfj9ctWnwzcysGDf4ZmYt4QbfrGKcKWZlqc1F23HrdRHV2Ts2bEWHJBjWBX/H9kvVPZnCDX5FNXF8ErNBeAa54fEpnRrx+CRmGWfvDGbWBl/SDZL2SNrasewKSbskbU5/53asu0zSdkkPS3pXWRU3K8qxbW3TTw9/DXB2l+XXRMSS9Hc7gKSTgGXA69JjPi/pkGFV1mzI1uDYthYZdBLzXpYCN0XEMxGxg2zs8NMK1K/ynFFRX47t5vJ+2V2Ri7aXSnofsBH4zxHxBNmkznd3bDPjRM/ACoA5x59YoBrj5bF3GsmxTb0z07xfdjdog38tcCXZdcQrgU8DH8xTQESsBlYDzJ1Y3Nf1yKaOp5OHs3dKN5bY7uQ4r74iXxzj3F8HytKJiMci4mBEPAdcxwuHtqVO9Owr8905e2d4xhXbnRznzTbO/XWgBl/SRMfdC4CpLId1wDJJcyUtAhYD9xarotnoOLatyQadxPxtkpaQfVntBD4EEBHbJN0CPAgcAC6JiIPlVN2sGMe2tc2gk5hfP8P2VwNXF6lUEwifaqk6x3b7tH2/9NAKJXGWgFn1tD25wQ1+Q3jmLGujoh2otmW4eSwdM2uttp3ecYNvZtYSbvDNzFrCDf6IeYwPMxsXX7QdMWfvmNm4uMFvOGfv1NMgHYC2ZZzUuZOUt+7D2l99SsesIdqWcWL5ucE3M2sJN/hmZi3hBr8inL1jZmXrZ7TMVwFrgePIThOujojPSjoGuBlYSDaq4EUR8YQkAZ8FzgWeBt4fEfeXU/3mcPbO6Dm2B9PtAmLb4jTPRdQqvTf9ZOkcIJvm7X5JRwKbJG0A3g/cERGrJK0EVgIfA84hGyt8MXA62QxCpw+rwotWrm9VJkJZnL0DVCy2zcrWzyTmu6d6MRHxFPBdsrk8lwI3ps1uBM5Pt5cCayNzN3D0tEklCnEmgg1L1WLbrGy5zuFLWgi8CbgHOC4idqdVj5IdFkO2w/yw42FdJ3uWtELSRkkbc9bZbOgc29YGfTf4kuYBXwE+EhF/37kuIoKcne+IWB0RkxExmedxZsPm2La26KvBl3QY2Q7xpYj4i7T4sanD2fR/T1o+ssme28DZO+VybA/H/Hlzxl0F60M/WToim/btuxHxmY5V64DlwKr0/7aO5ZdKuonsgta+jsNjy8nZO+VxbA/PxsvPGko5juty9ZOl80vAe4EHJG1Oy36XbGe4RdLFwA+Ai9K628nS1raTpa59YKg1tlJ12+EaPEaLY3tMFq1c7wSMHHrtl3n1M4n5t2Yo+8wu2wdwyQB1sYpq6o7p2B6fpsbUKA3yHvqXtmZmLeEG38ysJSoxHv4vnnAUe8ddiZrxz9ubr2W/erYRcA/fzKwlKtHDz8vj6Yxey7J3GmnUR4B1OULp9r7Upe551bKH7yv81eDPwaxeatngm5lZfm7wzcxaopbn8K27Xucdnb1jNj5V2i/dwzczawn38K0QZ+9YL3WeVS1v77sOrwncw7cSOHvHrJpmbfAlvUrSnZIelLRN0ofT8isk7ZK0Of2d2/GYyyRtl/SwpHeV+QLMBuXYtrYpMok5wDUR8anOjSWdBCwDXgf8AvANSf80Ig7O9CTCPcOy+L3taSSxPYimT3zjmCz+HpQ1PPJuYHe6/ZSkqYmee1kK3BQRzwA7JG0HTgO+PdPz7Fh1nrNJSuJJVLobVWznVZfzwUWUfY2nDrE9jPdAn8y3fZFJzCGb/WeLpBskvSIt62uiZ7MqcWxbG/SdpTN9omdJ1wJXkh2VXAl8GvhgjvJWACsAXv3qV+epM5B9g7ehJ1RXdcrQKDO25xx/Yu76OLbrpy5ZPQNPYh4Rj0XEwYh4DriO7NAW+pzoOSJWR8RkREwuWLCgyGswG1jZsV1u7c3y6SdLp+tEz5ImOja7ANiabq8DlkmaK2kRsBi4d3hVNhsOx7a1TZFJzN8jaQnZYe9O4EMAEbFN0i3Ag2RZEJeUkcVgxTlTwrFt7VJkEvPbZ3jM1cDVBeplI9D27J22xXbe88Z1joMqjV9TJf6lrZlZS3gsHRupOmXvtEXbe73jMK79wD18M7OWcINvZtYSlWrwmz5+SF34czBrpkqdw/d4OtXQ9uwdy8fXX+qjUg2+tVe3L5M2NySTV21g4+VnjbsahU1etYG9+58d+PHz581pxPvQrzydqvnz5uQuv1KndMwsU6SRrJKir6Mp70MZBnlv3OCbmbWEG3wzs5Zwg299c/aOWb35oq31zdk7VndtH2Nn1gZf0uHAXcDctP2tEfF7aXjYm4BjgU3AeyPiWUlzgbXAKcCPgV+PiJ1lVL4pmQzWXdk/P69ybMPs8b1o5fpCo522OQuqKeYcf+Ipebbv55TOM8DbI+KNwBLgbElnAJ8km+j5ROAJ4OK0/cXAE2n5NWm7UvgKvhVU2diG2eO75UNb2wBmbfAjsz/dPSz9BfB24Na0/Ebg/HR7abpPWn9mmmjCrFIc29Y2/U5xeEiaIGIPsAH4HvBkRBxIm3RO5vz8RM9p/T6yQ2OzynFsW5v01eCn+T2XkM3heRrw2qJPLGmFpI2SNj7++ONFi7MxqnMXt+zYLlxBsyHKlaUTEU9KuhN4M3C0pENTT6dzMuepiZ5/JOlQ4CiyC1zTy1oNrAaYnJz06cgaa0L2TlmxPXdiceVj2xdv26OfLJ0FwM/SDnEEcBbZxao7gQvJshmWA7elh6xL97+d1n8zIiof9FYf3b5IRO8vnl4c2/UwjI5D28bk6aWfHv4EcKOkQ8hOAd0SEV+T9CBwk6SrgL8Frk/bXw98UdJ24CfAshLqbfYiA7a6ju2WcEZfpp9JzLcAb+qy/Ptk5zynL/8p8K+GUjuzEjm2rW08tIKZWUtUrsEfZIxnq6Y6Z+/YcNRlf25LrFZuLJ2Nl59Vq+wO660J2TtWTN4LpeOKjbbEauUa/Lw8no5N6ZW9Ywbd46Nt2TuVO6WTl6++20ycM2kzaVv7UfsG38zM+uMG38ysJdzg28j5vPpw+H0sX9Pe49pftLX6aUtGRNnyDiVh+TUtVt3Dt8bLOyuQWVO5wTczawk3+GZmLTFrgy/pcEn3SvqOpG2Sfj8tXyNph6TN6W9JWi5Jn5O0XdIWSSeX/SLMBuHYtrbp56Lt1ETP+yUdBnxL0l+ldb8dEbdO2/4cYHH6Ox24Nv3v2/x5c1r3gwjLMiJG/EOp0mP70Jc1Lc+jXN73y9XP8MgBdJvouZelwNr0uLslHS1pIiJ291spj6fTTqPOiBhFbP/ziZezd2g1br66jL1TVwNNYh4R96RVV6dD22skzU3Lnp/oOemcBNqsUhzb1iYDTWIu6fXAZWQTPp8KHAN8LM8TexJzqwLHtrVJriydiHiSbL7PsyNid2SeAb7ACzMETU30PKVzEujOslZHxGRETC5YsGCw2psNiWPb2qCfLJ0Fko5Ot6cmen5I0kRaJuB8YGt6yDrgfSmj4QxgX57z92aj4ti2tikyifk3JS0gS67YDPy7tP3twLnAduBp4APDr7a1SYnZOyOJbWeelKdobNRlRq5hUZZwMOZKSE8BD3cuy/Nz+Gcf3b6px6r5MPQkCZdZkTL7jZED+/Zw8Ol9Y8mP7IztQYZ46BHblfocRlxmWeUOpcxRD+ORN7arMnjawxExOexCJW0cdrkus51lFjD02K7Le1bW51CXulaxTA+tYGbWEm7wzcxaoioN/uo6lCvpCmDHDOsnJX1ugKLLeP0vKlPSFZJ2dYwPs2qmB6fxZC6cqcwhqUuZg6rL61udYuSjvTYYIL5Hsl9L+nhHXB/suP2bg5Y5JJUrsxIXbesiNfj7I+JT465LXnnrLmkN8LUu48lYQ9U5vqdI2h8R83qsOzQiDoy6TlVSlR5+ZaXewyOSvgX8s7Ts1PSz+82S/lDS1rT8bZK+lm6/qLckaaukhenvodSDfkTSlyS9Q9L/lvR3kk7rePwXJX07Lf+3afmEpLvSc2+V9MsFX98nJN2Xylqdcs+nb7NK0oPpNX8qLVsg6SvpsfdJ+qUi9bDxaEF8/6mkayXdC/xXSVdJ+kjH+ockvTLdXq5s9NTNkj4vqXHtY+Ne0DBJOgVYBiwhy78+Na36AvCh9JP8gwMUfSLwabKf778W+A3gLcBHgd/t2O4NwNuBNwOfkPQLadu/Sc/9RrI8cSTd3HEo2/n3vo7yfqtj+bvSsj+KiFMj4vXAEcC7p70HxwIXAK+LiDcAV6VVnwWuiYhTgV8D/scA74ONUQPju5cJ4IyI+J1eGygbUuMC4F+k5z6U7L1plKqkZVbVLwNfjYinASStS8uPjIhvp9t/xrRGsg87IuKBVOY24I6ICEkPAAs7trstIv4B+AdJd5L9xP8+4AZlw/n+z4jYDBARv97H817T5XD9VyT9DvCPyMaN2Qb8Zcf6fcBPgetT7+5rafk7gJM6DgheLmleROzH6qJp8d3LlyPiuVm2eQfZF97GFNNH8OKB8hrBDX55DvDiI6jDO24/03H7uY77z/Hiz2T6BZaIiLskvRU4D1gj6TMRsVbSzaRD8mk+ExFru1VQ0uHA54HJiPihsnO4nfUkIg6kw/AzgQuBS8l6ZS8j6zX9tFvZ1niVj+8O/6+Pegu4ISL+yyxl1ZpP6czsLuB8SUdIOhL4l2n5U5KmJr7oddi3EzgZQNnMSIsGeP6lymZlOhZ4G3CfpH8MPBYR15GdRjkZsh5QRCzp8jfTzjAV7HslzSNr0F8kLT8qIm4HfovsMBvg68B/7NhuyQCvz8ar6fHdq96npHqfxguD4X0DuEjS/LTuWEmvHuA1VZp7+DOIiPtTz+I7ZOOl35dWXQxcJ+k54H+RnfZ4/mHp/1fIBtraBtwDPDJAFbaQjeA4H7gyIv6vpOXAb0v6GdnkHf2cw+wqIp6UdB3Z4GCP8sLr63QkcFs6GhDwn9Ly3wT+WNIWsji6ixfGnLEaaHp89/Bl4F+nC9F3A98HiIgHlE1x+Y10sfZnZPH8f4b8/GPltMwBdJ6rlrQSmIiID0v6NeBXI2L5EJ7jCmqeImf15PhuLvfwB3OepMvI3r8fAO+X9KvA1cAHx1ozs+Ic3w3lHr6ZWUv4oq2ZWUu4wTczawk3+GZmLeEG38ysJdzgm5m1hBt8M7OW+P9X1c8V61vVlwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a1 = BilinearForm(L2(mesh, order=order, dgjumps=False)); a1.Assemble()\n",
    "a2 = BilinearForm(L2(mesh, order=order, dgjumps=True)); a2.Assemble()\n",
    "A1 = sp.csr_matrix(a1.mat.CSR())\n",
    "A2 = sp.csr_matrix(a2.mat.CSR())\n",
    "fig = plt.figure(); ax1 = fig.add_subplot(121); ax2 = fig.add_subplot(122)\n",
    "ax1.set_xlabel(\"dgjumps=False\"); ax1.spy(A1,precision=-1)\n",
    "ax2.set_xlabel(\"dgjumps=True\"); ax2.spy(A2,precision=-1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Remark 3: Dof numbering of higher order FESpaces "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "* In `NGSolve` `FESpace`s typically have a numbering where the first block of dofs corresponds to a low order subspace (which is convenient for iterative solvers). \n",
    "* For L2 this means that the first dofs correspond to the constants on elements. \n",
    "\n",
    "* You can turn this behavior off for some spaces, e.g. for L2 by adding the flag `all_dofs_together`.\n",
    "\n",
    "We demonstrate this in the next comparison:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2oAAAEkCAYAAABe5sI+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xv4JVdd5/vPNze5RIhJNzEE2gaJYGRMgz8jTFC5iAMEgTliMDIYNGN7zqDi7UhQR6JP5jHoCKIi0gokMFyCYCQHGCTTmAEcBToQSCREMHQkMaRDTIR4IXbyPX9U/dL1q9619qradVmr9vv1PP30/u2qXXvVWqvW3rXru75l7i4AAAAAQDqOmLoAAAAAAICtOFEDAAAAgMRwogYAAAAAieFEDQAAAAASw4kaAAAAACSGEzUAAAAASMykJ2pm9jQzu87MPmdm501ZllWZ2X4zu9rMrjKzfVOXpw0ze72ZHTCzayrPHW9ml5vZZ8v/v27KMrbRsD/nm9lNZftcZWbPmLKMsczsoWb252b2aTP7azN7cfl8tu2TqqnGo1SOv1T6mpndx8w+amafLMvxq+XzDzOzj5Ttc4mZHTNkOSrlOdLMPmFm7564HId9xkzUT44zs3eY2WfM7Foze/wEfeSRlbH8KjP7spn99FzHRcamNMam8j2TGZ8Ymw4rxyzHpslO1MzsSEmvlvR0SadKOtvMTp2qPD15krvvcveNqQvS0kWSnlZ77jxJe939FEl7y79zcZEO3x9JemXZPrvc/b0jl6mrg5J+zt1PlfQ4SS8qj5Oc2yc5E49HFymN4y+VvvZVSU9299Mk7ZL0NDN7nKSXqziGHyHpdknnDlyOTS+WdG3l76nKIR3+GTNFP3mVpPe5+6MknaaibkYth7tftzmWS/o2Sf8s6dKxyzEGxiZJ6YxNUlrjE2PTVvMcm9x9kn+SHi/pzyp/v1TSS6cqTw/7s1/StqnLsUL5d0q6pvL3dZJOKh+fJOm6qcu44v6cL+nnpy5XD/v1LklPzb19Uvs39XiU4vGXQl+TdD9JH5f0HZK+JOmoRe014Ps/pPxQfbKkd0uyKcpRvtdhnzFjt42kB0r6vCSbshy19/5eSX8xdTkG3D/GpsPLNPnYVL7nZOMTY9NhZZjt2DRl6OPJkr5Q+fvG8rlcuaT3m9mVZrZ76sL04ER3v7l8/EVJJ05ZmJ78hJl9qgznyC4kxsx2SnqMpI9onu0zpdTGo0nbd+q+Vob0XCXpgKTLJf2tpDvc/WC5yljt89uSfkHSPeXfJ0xUDmnxZ8zYbfMwSbdKekMZcvVHZnb/CcpR9YOS3lo+nuO4yNhUMfXYVJYhhfGJsWmr2Y5NJBPpzxPc/bEqwhNeZGbfNXWB+uLFTwA+dTlW9BpJ36giVOFmSb81bXHaMbNjJb1T0k+7+5ery2bSPmgwdvum0Nfc/W4vQkceIul0SY8a+j3rzOyZkg64+5Vjv3eD4GfMSG1zlKTHSnqNuz9G0j+pFsIzZn8t5+A8S9If15cxLg5vHcem8r0mHZ8Ymxaa7dg05YnaTZIeWvn7IeVzWXL3m8r/D6iIRz192hKt7BYzO0mSyv8PTFyelbj7LeXgeo+kP1RG7WNmR6v4cHqzu/9J+fSs2icBqY1Hk7Rvan3N3e+Q9OcqwniOM7OjykVjtM8Zkp5lZvslvU1FiNGrJiiHpMbPmLHb5kZJN7r7R8q/36Hiy9FUfeTpkj7u7reUf89xXGRsUnpjkzTp+MTYdLjZjk1Tnqh9TNIpZZaaY1RcIrxswvJ0Zmb3N7Ov3XysIi71mvCrkneZpHPKx+eoiAnP1uYBUvqPyqR9zMwkvU7Ste7+isqiWbVPAlIbj0Zv31T6mpltN7Pjysf3VTEX5VoVX4ieO1Y53P2l7v4Qd9+poj98wN2fP3Y5pOBnzKht4+5flPQFM3tk+dRTJH167HJUnK1DoUWasBxDYmxKZGwqyzL5+MTYdLhZj01jTaprmGj3DEl/oyK+95emLMuK+/FwSZ8s//11bvtSdqabJf2bil8lzlUR77xX0mcl/S9Jx09dzhX3502Srpb0qfKAOWnqckbuyxNUXCL/lKSryn/PyLl9Uv031XiUyvGXSl+T9K2SPlGW4xpJv1I+/3BJH5X0ORXhJF8zYhs9UdK7pypH02fMRP1kl6R9Zfv8qaSvm6gc95d0m6QHVp6b5bjI2JTG2FSWJanxibFpS1lmOTZZuQEAAAAAQCJIJgIAAAAAieFEDQAAAAASw4kaAAAAACSGEzUAAAAASAwnagAAAACQmMlP1Mxs99Rl6Muc9kWa1/6wL4iRSt1SjsOlUhbKsRXlGEdK+5dKWSjHVpTjcKmUZZVyTH6iJimJSuzJnPZFmtf+sC9rxsyeZmbXmdnnzOy8yJelUreU43CplIVybEU5Wsp8bJLSKQvl2IpyHC6VsmR9ogYAvTKzIyW9WtLTJZ0q6WwzO3XaUgFYd4xNANoY9YbX27Zt8507d2557tZbb9X27dtHK8OQ5rQv0rz2h30ZzpVXXvkld0+nQJLM7PGSznf3/1D+/VJJcvdfX7T+tm3b/F++5ngdvKcYD486wvTNJz1grOIepqmNr735y6OWsUtfq5dRUi9lTqXfj1WOZfW47aivrlV9dC1HauNT27FJko66/wP9yAc8qHhc6wtb1htp3Eq9zXMvR/XYb7KoredaH6tIpSyLyhE7Nh21yhub2dMkvUrSkZL+yN0vDK2/c+dO7du3b5W3BJAYM7th6jIscLKkL1T+vlHSd1RXKGPGd0vSjh07dOzZr9yygX0XnjlwEdvbed57tvydQxnrUixzipbV437qMUqC49PSsUnaOj4d8/WP0Enn/HbUxjm+8rfs2N9EW+ctdmzqHPrI5XsAOXP3Pe6+4e4bKfziBgCbquPT1GUBMJ1V5qidLulz7n69u98l6W2Snt1PsQBgJTdJemjl74eUzzXaduwxCx+nJMcy5lDmFFGPs9V6bNoMd5QO7wtWWY9+MQ9N7Utbr6dVQh+jLt+HbFxwub50512Sik6375efukJxkDLaGiP7mKRTzOxhKr4E/aCkHwq9oNonNy64/N7wk3p/pS+HheqjWq8maXMWxrJ6rNZ5m9cBCWo9NgFYXyvNUYtRnwdStfnBW3+M+aGtMSZ3P2hmPyHpz1TMoX29u/917OtD/XXKvpz7cVQtszc83+frchbbB5GXLmNTNbFEve3X5XhYJ03jHW29nlY5UYu6fO/ueyTtkaSNjY3xUkwCWGvu/l5J7526HABQxdgEINYqJ2orX77fduwxW0KIMF+0NXJS76/V0Lsq07ihkLkfR9Xy10MYQ+GN1b/rqiGqkmYTlrqorXNue/Rjc57SouPBFjyH/ITGu+o6WA+dT9RWDS2SwnMZMC+0NXJS769N6ZJd44Yj5n4chcpfreMuIT5zCw/Mva0xjNAXeEKW5iGmHWnr9bHSHDUu3wMAAABA/1ZJzw8AayGUDpu06f2IreNQumraAnNEev71Qnp+VA2e9REAclOfdxYy1xC1IebeheahdXFC5vPQ1hW3W2jnm096gPZdeGbUusuOsbbH9KKxkLYDxsOJGgDUhOadrUuK5CHm3sWm2e+6DHlYx9stjCW2brvM/eRYHAfp+VFF6CMAAAAAJCaLK2pDh+BwyX541DdyUk+Nftuddy3MslVPK1/v2zn3+yFuBRBKz99Ux6YixHHR60zd0vOHwrlya6cchfoBVtOU2r2e1j82vXv9eKtuY9H2sTrS86MqixO1oUNwuIQ8POobOVmWnn9/Zb5IdVkovCi3fj/EyUpsev4qD7wutu7rQu2SWzvliBPh4TR9wa8/H5vevSn0LvZ90R7p+VFF6CMAAAAAJCaLE7UhUi6Txnlc1DdyFuq/XZdhqy4px+v1G1vfXV8HpK7vW4mEjhXSxQ+D9PyoyiL0cewQHPSP+kbOQv23PietPmdqCjnOjetSxtBrFrVFdV7UOiItPgDkJYsTNQDIQSpzn3KeG9eXrmnF51xfpMWfP9Lz54/0/KjKIvQRAAAAANYJV9QAYAXVcLKqrumw+9AmtX6OYZKbQmn2q2LSik/RTmPrOy0+tzlIz7I2ri7bDA0OrVe17DhqCq0lzLadpjbkVhbriRM1AFhBUwhK13TYfWjzRSjnMMnYsrdJKz7ntNd9f0FOJdQXh8S2cfXWFrEhdcuOI0L2+sGJLKoIfQQAAACAxMzuilpsGE/f6wFYT11CjYYeS9qEpLUJk0zNorLHtkWXkLDcx/9Q1kepuX6a9jtU/7n1pXVTbeP68+p52aL1AMSZ3YlabChM3+sBWE9dQo2GHkvahKTlfPIxRNm7hITlIjYjYOx+59x31l3TCVWbMOHYZV3WA1Ag9BEAAAAAEjO7E7VquEUo9KLv9QAgZMyxpP5ejGPxqvVjDc/nKrRvc95vHK5rX+iyjP4EdDe70MfYUIy+1wOAkOpYsnHB5Y3z1fqYF1t/zcYFl3cp8loK1XdojtcQr+vbsjICANIyuxM1AEhd7ByyvuZFMde2H7FzvPp63ZhyKCP608d8xdhl9Cegu9mFPgIAAABA7riiBgAjq6c1r4bGVZn6CYXMOQV/SkK3YgiFNzalQpe23rJBar6NwtCW3WYCeYi9BUOTVdLzn9DQh+hPmELTsZDbbVc4UQOQLTN7vaRnSjrg7o8unzte0iWSdkraL+ksd799qjIuUv9gqKaEr3L1E7aY+gdRLkL1uGpa/2W3URgafaRfU41Nq4awrpKenz6ElDQdC7mF4hL6CCBnF0l6Wu258yTtdfdTJO0t/waAMV0kxiYAK+JEDUC23P2Dkv6h9vSzJV1cPr5Y0nNGLVQHsSmvCRtKVx9py2nr+ZhqbOq7H3IbB+RqLreLyCL0sc0cjdh1+5j3AYyNfhvlRHe/uXz8RUknTlmYRertGDJmG9O/wpbN/2nrhBnXcSq3JEhM8mMTgLRkcaLWZo5G7Lqkq0aO6LftuLub2cKpFWa2W9JuSdqxY8eo5Qq145Tx8/SvsKFTms8J6f7DQmOT1H186tJH6cuYI+aoAUCabjGzkySp/P/AopXcfY+7b7j7xvbt20ctIIC1FDU2SYxPAApZXFFrk1o6dl3SVSNH9Nsol0k6R9KF5f/vmrY4h6u342133rUwo1o97fvQ4dz0r7BQCvumNlyWtrzv9PyphByS7n+hwcempltBhFLrx64XswxIRei2KNV1UpfFiVqbD5fYddc4Rh4Zo99uZWZvlfRESdvM7EZJL1PxJejtZnaupBsknTVdCRdblp5//4VnLlw2dDg3/SssNj1/VShteahtu7ZhKiGH696Xphqbmr6YLvvCGrNe12XAFGL6ZA79NosTNQBYxN3Pblj0lFELAgAVjE0A+sCJGgBo2pCxUMhhPYSsGipXLdfQYYtdQzAlNYb9zSkTYJdQv0Vt1tTW1WWh+so95DCV0M1cLWv/ZVdWQ+GNoeO5Hqbdx7G+LMsqGWoR0nQs5DYucqIGAJo2ZCz0JaO6LBQqN/QXlT4y6uaceWuZLvUfG0rZJiwy9y+sqYRu5iq2/UOhunXVUOzYbfZxrMf2BfoGFsl9LNxE1kcAAAAASAwnagCgrSEQ1vD81KplGbtcse9dX6+pXlOt41SE6nHO9ZXDcTgHoXru2tf6PtaHKCOQm9mFPvaRonrI7Q21TQCryeE4rJcxxfk8OdRjDkL1uHHB5feGmYXm7qTSJwAA3czuRK2PFNVDbm+obQJYP2PO52HcSkdsu+c4x4s5auMYYv5X0zaZowZ0R+gjAAAAACRmdlfU+k5RPUTK66HTaAOYj1CodDW0rcoWPLeq2NsEIF7oVgahel3W7mP1iSGE0sunGOqbq9g+1KbfNG1z0fZX2d4qZQRyM7sTtb4H6yEGfz5QAMSKTXcf8/wqYm8TgHihtg3Va5d2H6JPDCH2lgWERa4mtg+16Tcx6/axvVXKCOSG0EcAAAAASAwnagCQsFAa6qlSmZMaux9dU/DHpi2fW3r7Oe/b2EjPD+TB3Me7aLyxseH79u0b7f0ADM/MrnT3janLsYq+xqbQfLIxb8vBLUC6z/+a8zyoOe9bE8an9IT6obT49hKh9ebcf8cUGjPXZbwYU+zYtNIcNTPbL+krku6WdDD3wRAAVhE752joOTWk0u8+/2vO6eHnvG/IR5fbS5Cef3ix86Gp73H1kUzkSe7+pR62AwAAAADQDLM+TqVNqBFhSeOivjGW0K036qmmh0xvXy/HOoa8LWqLmNuihFKC516PXVLfE3KGvnW5vURoPdLz96N+3EvhdsI4Vj1Rc0nvNzOX9Fp331Nfwcx2S9otSTt27Fjx7dLVJtSIsKRxUd8YS+hLa9OXjyH6ZL0c65jWvOsJRCgleO6hg11S3xNyhr6telsR0vMPo+m4D62H4a2a9fEJ7v5YSU+X9CIz+676Cu6+x9033H1j+/btK74dAAAAAMzfSidq7n5T+f8BSZdKOr2PQuWoTapY0sqOi/pGCqZMLU5a83ikvg+nVmc8xaq6HGP0w+GFbhcyt/EuJ51DH83s/pKOcPevlI+/V9Kv9VayzLQJsyGmf1zUN1KQaj/Mfd5V39Zxn7u67c67tsy1lOJSq1PHABBnlTlqJ0q61Mw2t/MWd39fL6UCAPSmSzps5h+tl6Y+0rW/xG4D64f0/GkiPX+aOp+oufv1kk7rsSwAAAAAAJGeHwAGk8qtIapp2auWpWAe8hYCU6q3i0TK+abU/VUxKbublsXcCmBd634uYm/xQHr+NIVuabLoNh19Co3J6z5GcKIGAANJ5dYQ1Q+2pjTsm/ZfeOZh680t1CXULnPb11hNX36W9ZfQss2+FLvNda37uegSPquG9UJIzz+MKU+ACLtstmp6fgAAAABAz7iitkRs6FKbEKdUwqEwPNp6vS0KJVlm6D5TD0k6oRbytnm1IxTqknu/DoX4DB3Wk1sYT72/SIuvYCwKV2rqS01hlsvqPud6XAdNIY31ddSwHqGP6yt2nFnHtuZEbYnY0KU2IU6phENheLT1cMzsoZLeqCIDrUva4+6vMrPjJV0iaaek/ZLOcvfbpyhjly+PQ/eZekhSKCyyKXQt93499pf62Cx3KWoKYdvUJbyx6TXL5FSPOYxPfYsJQWwbPhuzHqGP+Vs2zsQsmytCHwHk6qCkn3P3UyU9TtKLzOxUSedJ2uvup0jaW/4NAGNifAKwMk7UlqjfqX3V9dqui7zR1sNx95vd/ePl469IulbSyZKeLenicrWLJT1nmhJ2M3SfCW1/iPEOW+vIGp5PVb2tu/SRvvY5p3qc6/gU0tQ+9bYKtWOXbaTY/mgnNM6kfqwPzdzHu5C4sbHh+/btG+39AAzPzK50942Jy7BT0gclPVrS37n7ceXzJun2zb+btBmbcp+fVde0PynNu2VuUpq63uag7/YMlePv3/CTB++65fqjW2+0R2OOT1OKTc8vhdO+t90G40D+1jE9f+x3J+aoAciamR0r6Z2Sftrdv1x89ym4u5vZwl+jzGy3pN2StGPHjuj3y31+Vl3T/qQ07zanuUnrpOttDvpuz1A57IgjJ/2eM/b4NKUu6fljl4XWYxzIH+n5mxH6CCBbZna0ii9Bb3b3PymfvsXMTiqXnyTpwKLXuvsed99w943t27ePU2AAa4PxCcCquKIGIEtl2NDrJF3r7q+oLLpM0jmSLiz/f1ef79sl5X5fuoaahTTtT/35UHjjotTKfeqazj1nsaFAU4YFdb3NQag9u4RFBlN7T5QmbqrxaUpN7Tp06GP1VhDrEjY3hDbhh03LutZ3aCxZlzG/CSdqAHJ1hqQXSLrazK4qn/tFFV+A3m5m50q6QdJZfb7plB/0XUPNQpr2p/58Nd16bIhSX9bxy1VsKNCUYUFd2yX0ulBa/ybB/jfdjZcmGZ+mNNVxGjs2rWPYXBtdww/7GI/WcYyPxYkagCy5+4fV/DXsKWOWBQCqGJ8A9IE5agCQia6p0od479hl6CY2XfXcUld3Scsdqiu/5+6DvRYQySG1ez/a1OO6jEcp4IpawuaWBnwqsfVIfSN1y1LkV+dp9JHmvLqN6uP6e6F/XcefVG5l0LUcsbeBGHpeJACkgBO1hM0tDfhUYuuR+kbO+ui/fRwrHDvTSuVWBkOUo0v696nT82N4pHbvx5Rz1NCM0EcAAAAASAy/NCVsyjTgcxJbj9Q3ctZH/62n125KeV1VT4dOGNq0UrmVQdcU/FJzWu7q301SSc+P8QRvz1BbD83a1mPssYjVcKKWMOZI9SO2Hqlv5KyP/lvdRijldVX9w5rvxdNKZRzrIwV/lzCqhNLzYySxtwdhbArrux6p734Q+ggAAAAAieGKGgAoj6yfY5ZxUShlbLhaU8hkqvWaonpbS/Opx2XhmW2XhdYjPf/8tRmrQplDu2RIjQ3jzeGYbVOPt915F6GPI+FEDQCUR9bPMcvY9UtFKGQy1XpN0Zwza475hdVe/sxPjvZmmESb/tQUdts1W2FsdtMcjtmu9diE0Md+EPoIAAAAAInhRA0AtDUrXqpZP3MsYw5lThH1CPSveuxUQ/OsYZ2u25v7Mdt3PaIZoY8zMMS8lRzm6/RtHfcZAAAAaeJEbQaGmLeSw3ydvq3jPuOQHNo/9zKmWuYUUY9A/5rmlDFHrZ2+6xHNCH0EAAAAgMRwRW0GFqVUTXGbqVvHfcYh9ZTh1RTz1TDYKUNk6+mkU7QsxfNmvbZJhx1KgT3XEOVFbZ162wOpqx5XoXVW3d7cj9m+6xHNOFGbgSG+qMz1y0/IOu4zDqm2f2yK+bFDO5pCTFISOo6aUmMvq8fY8KI5CbV1qm0PpC7m2GlzfDWtO/djtu96RDNCHwEAAAAgMZyoAUBNKK3ylCmXc0/3HEplPcTrckZ6fqB/pOfvB+n5x0PoIwAAAGYvdopDaF6stHXebYzb7rxry7zn+jZym3ebQxnnghM1AKiJTY0+9hypHNLzh3Sda7aOc9RIzw9MJ3bMiU1N33UZQOgjAAAAACSGK2prpE1a8dh1p0xVDgylnmK+2s+rTOMeA7nfQqJ+C4RquE8o1CiUCropnCj38Sh0m4Mc2x7ISWza/dDrm9aLWQZs4kRtjbQJm4pdN/dQLGCR+hf8alr5Kte4x0DOJx7SMKn7m9bLfTzKva2BnMWm3W/7+lWWYT0R+ggAAAAAieFEbY20SRUbu+6c088Cm9Y1BfOYYus4lAqatgDQh77HozbLgCpCH9dIm1Ca2HUJz8GUzOw+kj4o6WtUjGfvcPeXmdnDJL1N0gmSrpT0AnfPOxYOQDYYmwD0gStqAHL2VUlPdvfTJO2S9DQze5ykl0t6pbs/QtLtks5d5U1CqZqZp9mP2DoOpcOmLZCQUcYmDKPv8ajNMqCKEzUA2fLCneWfR5f/XNKTJb2jfP5iSc+ZoHgA1hRjE4A+EPoIIGtmdqSKEKJHSHq1pL+VdIe7HyxXuVHSyW22WU+5H0rVfELmKfObDHHbgS4p+GPTYZPCPh+hfjCncPohxiaspqnv9T0ehdarj1X190az+ueS1FyPcxlLlp6omdnrJT1T0gF3f3T53PGSLpG0U9J+SWe5++3DFRMAFnP3uyXtMrPjJF0q6VExrzOz3ZJ2S9KOHTu2LIsNoXPNd57mEGGEXcJ9ln0h2n/hmasVCqNbl7CvrmOTFB6f0F1s2GKTVdLzM1atLvS5NNexJCb08SJJT6s9d56kve5+iqS95d8AMBl3v0PSn0t6vKTjzGzzh6iHSLppwfp73H3D3Te2b98+YkkBrJO2Y1P5GsYnAMuvqLn7B81sZ+3pZ0t6Yvn4YklXSHpJj+VCC7EhSm1CmYYIe0Iz6rsbM9su6d/c/Q4zu6+kp6qYrP/nkp6rIrvaOZLe1Wq72hoSI4XDYOZoUR30uc3681phGfKyrB/MwVBjE1bT1Pfq66hhva6hj3Pq21Nax8/mrnPUTnT3m8vHX5R0YtOKXL4fXmyIUptQJrKnjYv67uwkSReXc0GOkPR2d3+3mX1a0tvM7AJJn5D0ujYbbQqJCa03N7F10HWbMc+vsgzp6tIPMjTI2ITVxPSxPsac+noz69uTWcfP5pWTibi7m1ljnbj7Hkl7JGljY2NOdQdgYu7+KUmPWfD89ZJOH79EAMDYBKAfXdPz32JmJ0lS+f+B/oqEtqpZgkIZg2LXa7suVkd9p6XeHtW/rWG9uRmiT4bqsesy5If2xFSa+l4f41FoPfp2P9bxs7nrFbXLVMRWXyhirCcXO5+pzbwn5kiNi/pOS709Ni64/N7HJ9TmEIbmF+Y893CIsi6bF9v2dRsXXK6d571H0uL0/LHzdWNfh35QvwAQJyY9/1tVJA7ZZmY3SnqZihO0t5vZuZJukHTWkIUEgCmF5hB2XYatuqRsj637vl4HIG9d0vPHLgutx7jSj3VMzx+T9fHshkVP6bksAAAAAAD1kEwEAOYulKq+6zJs1SVleyhV8xCvA5C3bccec+/Vlvo44JV1brvzrk63C6luQ9KW92oK066/ru/Q4FCo99Dv3bdq+y3bl7ngRA0AlgilBO66DFutmrq/TTrsrq8DkLfYE5HNk6q6RePD/gvPbLW9sUP25hQumPqJ5BC6Zn0EAAAAAAyEEzUAWCKUZrnrMmzVJc1yKFXzEK8DsB5i0/PHjhdTppVfx5T2c0Lo4xoZIlX4lOnHc059jryE+lY9HX99HsIUcjw2upRxWbr/0JyQdVTtF7nNTQGAdcSJ2hoZIlX4lOnHSX2O1KSS9p1jo/u8jDnXV5dbIADrKPZYiT12ppwnNqc5auuI0EcAAAAASAwnamtkiHkYU87tYF4JUmO1x/W/pyrHOgq1RX29damvUB0AOKTv8aLNeNS3Kd8bqyP0cY0MMQdhynkNzKlAalJJ+85tAeLrYJ3S83e5BQKwjmKPldhjp+t41Icp3xur44oaAAAAACSGEzUA6Ekqad8JC26Xknpd6ou03EAc0vMjFYQ+AsAK6inPU9AmLDjHVP6b6mWXSMEf0nfbhuo/t76Usmtv/nLwVhPV2yzELqN9wvqum65jcrUNu95Sg3YOC922RJr+mOFEDQBW0Hca57HlnMq/a9rpHNsbNW4FAAAgAElEQVQpRancjmLuDt5zqMfSz+ev6TOFVPrDSP0znNBHAAAAAEjMpFfUYkNucg7NQYE2HBf1PZ5qqET9eWlr2EqK6qE1OVlU9pi2CC3LrQ6mFKp/6nEY9PP5a/pMqa+DfqT+GT7piVpsyE3OoTko0Ibjor7H03ca57HlnMq/j7TTubRTilK5HcU6oZ/PX0zb0H79Sf0znNBHAAAAAEjMpFfUth17zGEZo1ZZD+miDcdFfY+nWtfLMrBVM7cNGY7aJhtfzn1lUdlXzYZXbaeuWdZSFZvdLHa/Q/WfW1/KBaGP80fo47gIfQwgtej6oA3HRX2PJ7auN7/8S8OHo7bJxpdzXxmi7NV2mluWtdjsZrH7nXPfyRWhj/NH6OO4CH0EAAAAALTCiRoAjKAaCjZ0WFj9vcZ879xV68cans9VaN/mvN+5O+qIQy2yrK1il9Gu6WpqQ47LYaR+zHDDawAAgER980kP0L4Lz5y6GBhJbEhx1zmnoWWEM6eHEzUAWTOzIyXtk3STuz/TzB4m6W2STpB0paQXuPvkE4zGvGVCmzlq2Cp2HleO+p6jhuVyGZ+Qnz6OZ471+Hqcqn4IfQSQuxdLurby98slvdLdHyHpdknnTlIqAGB8ArACrqghKfW04qHL8G3WXTfrUjdm9hBJZ0r6b5J+1sxM0pMl/VC5ysWSzpf0mkkKWFFP49+Uqr+PtqunTb/tzrsmTzGci9DtFkKhRrFj1ZShRqF9q/aRKvpLdzmNT8hPbFr5tsvWDen5gRbahIeNGUqWmzWqm9+W9AuSvrb8+wRJd7j7wfLvGyWdvOiFZrZb0m5J2rFjx8DF3DrvIJSqv4+2q3/5b0o5j8OFTpy6pu5PJZwydt+q6C8ryWZ8Qn5i08r3sWzOSM8PAAMws2dKOuDuV3Z5vbvvcfcNd9/Yvn17z6UDsM4YnwD0gRM1JKVNKtQU0qamak3q5gxJzzKz/Som5z9Z0qskHWdmm9ECD5F00zTFaxZqnyHabk36w+C6prDPIfV9DmXMTLbjE/LQx20aONZJzw+00ma+xlznXfVhHerG3V8q6aWSZGZPlPTz7v58M/tjSc9V8eXoHEnvmqyQANYS4xOAPnBFDcDcvETFxP3PqZgT8rqJy3OY2PT5fc1hWqM5i4PqOtcslTlqITmUcSaSH5+Qh9Ax23XZOoqtq6nqhytqALLn7ldIuqJ8fL2k06csDwBsYnwC0BUnagAwsnr6/Gr69irTMOn60U3X1P1N6Z+lrbdpkLRyW3cV2jcAq6mP45Iaj7fY4570/P1YNvY13Vanvmyo265wogYAIwulz69yDZOuH90Mkbq/ab2UUvcDWE1oHO8afkh6/n7Ejn2h2+oMGULKHDUAAAAASAxX1LCSPsKyum6v7/deV9Tj9GJDWNY1NCUHXcKQQuvR1sB8VMeHvsIPCX0c1xBtGIMTNayk74w4bbaXQjaeOaAepxcbwrKuoSk56BKGFFqPtgbmwxseh9Zrs83YbRD62N0QbRiD0EcAAAAASAwnalhJ33dtb7O9FO4YPwfU4/Sq9W6152mfPMS2YXUZbQ2sh/qxHRovumwzdszp673X0RBtGIPQR6yk7/lMbbbHXKp+UI/jW5SqucmY7cN8xbBQCv4uTphxHYfqaq77DAB940QNAEY2RKrmPjBfMaxaJ6F26rpsTmLrClgHQ4z5jEfjmupzm9BHAAAAAEgMV9QAYGRTpfldZlG5cAjpsOMtqyusr2UhxG3DrxeFkqcWdrvt2GOiy9hlm6F6jF1mOnRT51TrcUpDtGGMpSdqZvZ6Sc+UdMDdH10+d76kH5N0a7naL7r7e3stGQDM1FRpfpeJLde6Ih12vFVvV4D5ig3Ziw0hSzWUvGqIE5u+t7l5kialW49TmurkNCb08SJJT1vw/CvdfVf5j5M0AAAAAOjJ0hM1d/+gpH8YoSwAsBamSvO7DOnhw0iHHW8d9xlx+r4dSarjaW6oxzStMkftJ8zshyXtk/Rz7n57T2XK0hBprUmV3Q/qEamp98GNCy6/93E9ZXts/+2jn3NshIXqp9qGsa/buODyxjkhXdswlbT49CUAWF3XrI+vkfSNknZJulnSbzWtaGa7zWyfme279dZbm1bL3hBprUmV3Q/qEakL9dHY/ks/n1aXdPT1NuujDUmLj9SF+miXYyB0HHEMxKMe09TpRM3db3H3u939Hkl/KOn0wLp73H3D3Te2b9/etZwAAAAAsDY6naiZ2UmVP/+jpGv6KU6+hpjbwXyRflCPSJ01PF62rMt6GEZTnbdpsz7asEs5gDGF+miXYyB0HMW8LwrUY5pi0vO/VdITJW0zsxslvUzSE81sl4qrofsl/fiAZcxCDqlX1xX1iNSF0uJ3SeVPCvTxrZq6P9TuQ5cDGFNsH43ts6ne7iQ31GOalp6oufvZC55+3QBlAQAAAACoezIRAEBPQuG5saG7hPhOq0sq61A67K5tSEptpI70/GmiHtNk7uNdxNzY2PB9+/aN9n4AhmdmV7r7xtTlWMXGxobrOb+eRFrzkClvNdH1NgGSFtZrqnWcilA9zrm++r69wFzGJ747NYsdcxiP+tN0nMbW9zqNaU1ix6ZV7qMGALORQyriKVPw93GbgKY5ECnVcSpib9kwNzkch0hL7JjDeNSfpuM0tr7XaUxbFaGPAAAAAJAYrqglrGuo0bpdPu4L9b3eqmEZ9edTUQ8xSfG9F623LMA+pTpORage51xfORyHSEvsmMN41J+m47S+jhrWW6cxbVWcqCWsj1AjxKO+82Rm+yV9RdLdkg66+4aZHS/pEkk7VdxC5Cx3vz20nRzSmk+Zgn/I2wSkVMepGCJ1fw5yOA5j9TU2IazvtPI59rWxxdRRm7agzpsR+ghgDp7k7rsqE3PPk7TX3U+RtLf8GwDGxtgEoDOuqCVs27HHHJbJaJX1EEZ9z8qzJT2xfHyxpCskvST0gmq7LspOlYJ63+s7Q15IbOjjouMjNjsYDllWjzvPe8+9j2OzrOUQpr0GoY+txyaEEfo4PkIfx8OJWsJiP1Rz+PDNAfWdLZf0fjNzSa919z2STnT3m8vlX5R0Yv1FZrZb0m5J2rFjh27IoF3rfW/zy7o0fNay2PAijo9+hOoxtt1zzGQ3p9BHdRybpMPHJzQj9HF8hD6OhxM1ALl7grvfZGYPknS5mX2mutDdvfyipNrzeyTtkcr7qAFAvzqNTeUyxicAzFEDkDd3v6n8/4CkSyWdLukWMztJksr/D0xXwuFUQwat4fkh3otQxWmF2n3MPjGE3Mtftc5j05jqY1OX4yP3vja22HqMbQvqvBlX1ABky8zuL+kId/9K+fh7Jf2apMsknSPpwvL/d01XSgDrhrFpPIRb5++2O++6N6S77dzmnOfkxuBEDUDOTpR0qZlJxXj2Fnd/n5l9TNLbzexcSTdIOmvCMg6mOu9o6PlI3JYiHbHtnuMctTH79MDWemzCvDUdp13Ho9hlcxvvYnCiBiBb7n69pNMWPH+bpKeMXyIAYGwC0A9O1IBSNdV56LJ57HpDbRPrJdQ3xkxlXr99QTVMhf7aTb1tJUWNA8vaPef09qHbZIx5OwoAzfpOzx+7rO025oATNaAUG9rVJgRsiG1ivYT6xpipzKtfhKvp4emv3YXaNlSvXdo9l7SBfdyWAMCw+k7P33bbfW0jB2R9BAAAAIDEcKIGlGJTxbZJKTvENrFeQn1jqlTm9Nd+hNKKx44Xc0vPHzLnfQNy0nd6/q63UViHMYHQR6AUO8ehzVyIIbYJAACA+eNEDQASFjuPacw5O8yp7EfXOWpzTs8fMqPU/UDWSM8/HkIfAQAAACAxXFEDgJ6EUul3vQVDNV35ojlqi9KVD50+P1SmVHRNfT9mCvhF9RhTr6EU9tVtzC29/bL9xryE+qsUdzxjGKHPnmo73XbnXYOl5+863uU29nGiBgA9iQ1laxOWEfoQaVo2dPr81D/YpOHDCvvQtR67vi739PY59Dv0J/ZYzKX/zknssVgdc6pWSa2//8IzW793bPhkigh9BAAAAIDEcEVtAl1DoJCGNu3Xd1vTd9JWD6+IXTZmOdbFojqIqZPq6+rP527O+4b5WdZf132My0FsG7Zdtsp7d93eVDhRmwAZ0/LWpv36bmv6TtqawiuWLRuzHOsiti1Cr4t9TS7mvG+Yn9j+Sv9NV5cxZ5WwyLbr5tB3CH0EAAAAgMRwojaBaoYqslXlp0379d3W9J20hdpnzLajnxxeB7F1Ul1mDc/nas77hvkJ9VfGuDzEtmGbZau+d25jH6GPE2BeUd7atF/fbU3fSVts+9x2512Dps+vby/3tOxdDLFfuddj7HzaUCpribToADAWTtQAYARjpn1P6b1zE1tXc6vHpv0mLTqmQnr+/PUxnnYda7uMaSki9BEAAAAAEsMVNQAYwZSp0UnLHm+VdNI5a5PKmrToGMO2Y49pDDWWtobgTiUUDi0RJtxlPK3XXT0UuzptoGm9zXluseHcMeq3Rwq9d59tzYkaAIxgytTopGWP13c66Vx0SWU9h/1GunI4sSE8M6zNeLr/wjOXbm/zJE1aHiIZs702QrdHGjKcktBHAAAAAEgMJ2oAMIIpU6OTlj3e0CmjUxWbypq06MAh3EIgrO/6Cd12Zegxear3JvQRC9VjcXMIQRjLEHUzZX3T1uNItV5zTznft3Xc567qt5mQmueVrHu/AoAuOFHDQqFY3HU3RN1MWd+09fyta8p5xOuSyjp2Gf0Kc8UctbC+62eqeWJTvjehjwAAAACQGE7UsBCx1c2GqJsp65u2nr+mNOq2ZBnWR0x7L+svsdsA5iL2eFjXft93/dRfM+bn11TvTegjFmIOQbMh6mbK+s65rc3sOEl/JOnRKqIPflTSdZIukbRT0n5JZ7n77RMVMQnrmnIe8bqk52/7+jbrzQHj0/zFjq3r1O+r+q6fppDq2PddxVTvzRU1ADl7laT3ufujJJ0m6VpJ50na6+6nSNpb/g0AY2N8ArASTtQAZMnMHijpuyS9TpLc/S53v0PSsyVdXK52saTnTFPCdKxrynnE65KeP3bZOvYrxqf1QHr+MNLzr87cwxfpzOyhkt4o6UQVV/T2uPurzOx4tbx8v7Gx4fv27euh2ABSYWZXuvvGBO+7S9IeSZ9W8Wv1lZJeLOkmdz+uXMck3b75d5OhxqYUb31QL5OkhWVMsex9ia0DhOuqnna/uixUj2PeEoLxaTVt2p/jBl0N0c+axpmh+23svvz9G37y4F23XH/0su3FXFE7KOnn3P1USY+T9CIzO1VcvgcwraMkPVbSa9z9MZL+SbVxyItfohb+GmVmu81sn5ntu/XWWwcpYIq3PqiXqamMKZa9L7F1gHBd1VNSx9ZjbMruzCU/PsVo0/5AV0P0sy63HOlD7L7YEUdG5QlZeqLm7je7+8fLx19REWN9srh8D2BaN0q60d0/Uv79DhVfjG4xs5Mkqfz/wKIXu/sed99w943t27ePUmAAa4PxCcDKWmV9NLOdkh4j6SOSTnT3m8tFX1QRGpmNOYf1IB1997Mp+21qx4y7f9HMvmBmj3T36yQ9RUWY0aclnSPpwvL/d01Vxnq4RQoWlWlRGVMse19i6wDhuqqvpwXrLtvmom3MQQ7jU4y27Q90MUQ/axpnum6vy/uG9iU2PWT0iZqZHSvpnZJ+2t2/XIRWl+/l7mbWePle0m5J2rFjR+zbDY5QF4yh7342Zb9N9Jj5SUlvNrNjJF0v6UdURAq83czOlXSDpLOmKlxsOt8xhcrUtCyVsvcltg7QPSV1l/TVM6z7pMenGFOmQ8f6GKKfrXrLka6iPzsjzxKjTtTM7GgVJ2lvdvc/KZ++xcxOcvebl12+VzGhVhsbGxzLAHrj7ldJWpQo4CljlwUAqhifAKxq6YlamZXodZKudfdXVBZdpowu39dtO/aYw7KyAH3ru59N2W85ZtrrUmdDh5jWwzJOqJTRJO087z2HrVcve2phsG0tapch+/aYWQ77FhvGsyi7WVNf2lbrc039rC7neswVoY8YA6GPzWKuqJ0h6QWSrjazq8rnflHFCVo2l+/rGNQxhr772ZT9lmOmvS51NnSIaT0so1rGzS/Wm8v2X3jmwm0kGgYbbey+nHOmvGVhPE19pN6XNn3pzrsaX7NMzvWYK0IfMQZCH5stPVFz9w8HNsflewAAAADoWcx91AAAI6mGfw0Rhhfafux7D13GuanWkTU8n6p6W3fpI33tc871mKtQ+9MG6MsQ/axpG0P329h98XvuPhizvVbp+QEAALAeCHkfX9NczEXzOUNzPXOaRzxE+WK3GZr7Kk1fj5yoAUBChp7/Fdp+7HvnPkdtbDnPrYrtL6HX9bXPOdcjEKupn4f6fGgZx0dY7LgSW4+hMbO6fTviyKhzMEIfAQAAACAxXFHDKGLTeeee9jsVQ9QjbVMYuh6Gvg1CffuhMJvNzH31/QxtY4g6yT0te9d09CkI3cqg3kfqy5pCtrq2Z1M9hspx9IMedlrbfQam1CatfCiFfT1FPBZrqu+u9ThFen5gZYRUjWuIeqRtCkPXw9AnHvXtN6Xkrz5f38/QNoYO18wx5C2Hk8kmobKH2j10e4dQ6v4uZQmVIza8CEjFqmnl68u4dUJYU/10rce+0/MT+ggAAAAAieFEDaMg7fe4hqhH2qYwt3po2p82+znmLQVICZ6OVFL3h8oRmwIbSEVsWvnQcTS3z6kh9V2PpOdHlmJDf3IOEUrJlKluAQAAsDpO1ACghbnN1Wvanzb7OeYtBXKcozZXqaTuD5WDOWrIDen5x0V6fgAAAABAK/zSBAAtDJ0+P6SeBl/Slr+7hKc27U+bFPyL0hH3Kef09l2F2jqU+n7MEOVQ6v5lc9Sa2rNL6v5gOmxS3iEzodtQ1Me+0LjY9jYaOdzuJHZcbLMvyz5f2tZjFen5AWBkU36QdQ01C2nanzYp+KPTEXeU+peHIcSGz0wZCtq1XWJT/sfuW7D/cRMpZKbv8S527M4hRDJ2XGyzL7H1HVuPMc9LIj0/AAAAAOSKEzUAyETXdOhDvHfsMnQTm+J5brcr6JK6n/T8QLPYsSSH8WPKfek6JpOeH2srNGcmN1PuS+x7z6m+AQAAUseJGrI1pzTpU+5L7HvPqb5zNcQctb7fm77RjxzmqA2hS+p+0vMDzYaY1zWVKfeljzGZ9PwAAAAAMAP80oRsTZkmvW9T7kvse8+pvnMVSodeTx3cJTQ1FN5aT2Fcfa+h0/Ovo2Vt3ZSyu0t6+yF0LUeX1P1VpOcHtgrevqK2Xuqm3JfY915UDtLzYy3NaY7UlPsS+95zqu9cxaY17xr2EQrtqL53lzTFaKfr8dYlvf0QuoQwSt1S91eRnh/YKvb2KTmM3VPuS5fb0JCeHwAAAABmiCtqADADfYQfxm4jFALCBYxpLQsJzLkcTdtctH1CH4FCbBh1DlMa2uxLU6j0on0O1ceiKABCHwEggpk9UtIllaceLulXJL2xfH6npP2SznL328cu35i6hGV03UZoPb4XTys6JHBgQ5Qj5rWphD4yNiEVc5q20GZfmkKlY7MyhsK0CX0EgAjufp2773L3XZK+TdI/S7pU0nmS9rr7KZL2ln8DwCgYmwD0gRM1AHPxFEl/6+43SHq2pIvL5y+W9JzJSjWSathK1xCW2G3U1+vjvdGPav1bw/O5lqNpm/XtV9fze+4+2PkN+7PWYxMwhS7jRR/LYteLHZsIfQQyEkqf3mW9obY5kR+U9Nby8YnufnP5+IuSTpymSOPpoz2q29i44PItKfilrXH8SFMqx2XsmNN2vkim1npsAtAdJ2pARkLp07usN9Q2x2Zmx0h6lqSX1pe5u5vZ4VNXzHZL2i1JO3bsGLyMuQm1dyiOP7W+gfTEpu6PXRZaz444ctLvOV3GpvJ1jE/ACprGmT7GnK7jUXVZ7NhE6COAOXi6pI+7+y3l37eY2UmSVP5/oP4Cd9/j7hvuvrF9+/YRiwpgjbQemyTGJwAFTtSAjHSdQzT2Nidwtg6FFknSZZLOKR+fI+ldo5coc/WY+6Zwx/oywiKxTGxfil0WXG/6NKSMTcAEOo0XPSyLXo/0/MD8xM4/aTNPZYhtjsnM7i/pqZJ+vPL0hZLebmbnSrpB0llTlC1nXVL1L1sXkLql7u+cDnvCXw4Ym4DpdLqdRw/LoteLHJs4UQOQNXf/J0kn1J67TUWmNQCYBGMTgFWZ+3i/f5rZrSp+QaraJulLoxViWHPaF2le+8O+DOcb3D3rSRSVsSmVup28HEc/6GGnbU523kwjfO/kZ9e9vwbWl/k9dx/8twOf/+QARZq8TkqUY6vW5aj2rVBfil0WWu/gHV/0u//lK1lP8yjHp39SGu0tZdz3BkI5tkqiHE3jTB9jTtfxqLosdmwa9URtYQHM9rn7xqSF6Mmc9kWa1/6wL4iRSt1SjsOlUhbKQTmmkNL+pVIWykE5lkmlLKuUI+tfmQAAAABgjjhRAwAAAIDEpHCitmfqAvRoTvsijbw/VviAmT2gh+38jpl9zsw+ZWaPlbTHzLab2ft6Ku6U5tbPUpJK3SZdjh6P1UeZ2V+a2VfN7OdjymJmV5jZlKEsm+W4c6w3NLMXmtnvNZTjvWZ23JLXD1lnN0e0Xb08zzSzX+u5HKkcM0NJaf9SKQvl2IpyHC6VsnQux+Qnau6eSiWubE77IrXfHzNbNYvoMyR90t2/vOJ2ni7plPLfbkmvKW8eequKLxVnrLj9Sc2tn6Uklboduhyxx2qgHH0dq/8g6ack/fcVyjK4an11KUd5Ytv4eWtm+9tuc7Mc7v4Md7+j7eu7aOg3V3bYxnskfZ+Z3a+XgimdY3coKe1fKmXpsxyr/PhULUfLH5+WlekiM3tu+Xjpjy2b5Wj4Yae+7e1m9hEz+4SZfWeLMr3QzB4cU442zGyXmT2j8vf5PdTf2yT9+Srb6MsqfXXyEzV0Y2Y7zexaM/tDM/trM3u/md23XLbLzP6qvJp0qZl9Xfn8FWb2cjP7qJn9zaKD08webGZXVf7dbWbfUB7U7zSzj5X/zijXP9/M3mRmfyHpTWZ2HzN7g5ldXQ4ATyrX+5byfa8qy3XKgt16vsqbf5b79xkze3O5n+9o8aH+bElv9MJfSTrOzE4ql/1p+T7AKDhWm7n7AXf/mKR/61i3Z5flv8bMXl4+9wNm9ory8YvN7Pry8cPLfa9vI9QGv21m+yS92MweZsUXsKvN7ILaNv7fsq4/ZWa/WqmX68zsjZKukfTQLvtY8WAze5+ZfdbMfqPy3vvNbFv5+L+W7/lhM3urbf2i8wNL+pOZ2W+WdXm1mT2vfP6JZvYhM7tM0qfL536p3M6HJT2yso1vLMt4ZfmaR5XPX2Rmf2BmH5H0G15kMbtC0jNXrBNAUlI/FEf/+DSxp0i62t0f4+4favG6F0oKnqh1tEtFG/TCzI6U9BpJv9DXNifj7vzL8J+knZIOStpV/v12Sf+pfPwpSd9dPv41Sb9dPr5C0m+Vj58h6X8teY8XSXp7+fgtkp5QPt4h6dry8fkqflG9b/n3z0l6ffn4UZL+TtJ9JP2upOeXzx+zuX7t/W6Q9LWV/XNJZ5R/v17Sz5ePXynpqgX/ziuXv3uzrOXfeyVtlI9PVjE4Td6G/FuPfxyrzcdqZXvnb74moj6vkLSh4svC30naruKeoB+Q9BxJXy/pY+W675D0sfK4P0fSry/YXqgNfr+y3mWSfrhS33eWj79XRViLqfjx892Svqusl3skPS5in/YvWf5CSddLemDZRjdIeujma1Wkw/72sm7vI+lrJX220g5L+5Ok75d0uaQjJZ1Y1u1Jkp6oIjX8w8r1vk3S1ZLuJ+kBkj5XeZ+9kk4pH3+HpA+Ujy8q6+XIyvs9X9LvTn188q/9v7JvXyvpDyX9taT369C4skvSX5XH1aWSvq7SB18u6aOS/kbSdy7Y7oO1dZy4W9I3lMf4O8tj+WM6NNacL+lNkv5C0lvLvv+Gsn9+QtKTyvW+pXzfq8pynbLgvd8i6YmV/fuMpDeX+/kOSfdrWUfnK35M+5Vyv65ROZaUz18k6bmV+tsIbONHynr9aNkuv1fZlw+U+71XxWfCrvL4vrWsk/uX73VNWXc/0/Aez5V0p6TrytfdV8UJ3yfK171e0teU6z6jrMMrJf2OpHeXz9+/XO+j5euereJzplqe55X19/pyv6+X9FOVcvynSnu+VuW4UpbttyR9UtITVIzHn5d01NTHzCr/uKKWt8+7+1Xl4ysl7TSzB0o6zt3/d/n8xSq+NGz6k+r6TRsuf4X/MUk/Wj71PZJ+z8yuUvGF5QFmdmy57DJ3/5fy8RMk/Q9JcvfPqPhC8U2S/lLSL5rZS1Tcd2tz/arj3f0rlb+/4O6bv4D/j3LbcvefcfddC/5d2LQ/FQc0zK9BQAjHavtjdZlvl3SFu9/q7gdVfKn6Lnf/oqRjzexrVVzFeouKev1OSVt+OY5og0sqj89Q8WVQKr4cbvre8t8nJH1cxUnv5lXIG7y4qn8YM3v15tVQFVfLNq+M/lLD/u519390939VcWXrG2rLz5D0Lnf/17Jt/r/a8mX96QmS3urud7v7LZL+t4o6lqSPuvvny8ffKelSd/9nL64+XFbuz7GS/r2kPy736bUqTvQ2/bG73135m7E4b6dIerW7f4ukO1Sc6EvSGyW9xN2/VcWX95dVXnOUu58u6adrz0uS3P3vN8cIFScb73T3GyS9StIr3f3by/f5o8rLTpX0Pe5+toofUNzd/52ksyVdbGb3kfR/S3pVud0NSTcu2J8ztDWM95Eqfqj5ZklflvRfJMnMXlmLZNj8d15ctS30e+7+7e7+aBUnP62uNFsRMfSr5T48QUWdbPpdSReX7fFmSb9Tfhb9iqRLyjp5lKST3f3RZd29YdH7uPs7JO1T8UPeLi56rUIAAAbMSURBVBU/0F0k6Xnl646S9P+Udf5aSU93929TcaK96ZdU/IBzuqQnSfpNSUdXy+Pum+PuoyT9B0mnS3qZmR1tZt+s4kTujLIMd+tQlNT9JX3E3U9z9w+7+z0qfkg6rU19pmbVS8WY1lcrj+9WcYDHvuZule1vZm+Q9BhJf+/uzygP+tdJepa7b06YP0LFL8P/Wt2YmUnFr61B7v6WMuzlTEnvNbMfd/cP1FY7aGZHlAeXVAwCWzZTvucrVRzgdW8rvwDepK1hRg8pn5OKX9wWffEEhsSxutXmsTqU/6PiF+brVJyc/aikx6u4ithGvb4W3XjUVFype+2WJ812Lnj9oQ25v6iy7v7yS0dIvQ+1/fw+rD+1sLTfqOh3dwT2o74NxuK8xf749MeV17T98ekJ5VPfI+nUcgyTwj8+/a5U/PhkZtUfn37JzB4i6U/c/bML3nbZj08/Jem/u/vPNJV7BU8ys19QcZX6eBVXKes/tIR8h8ofrSTJzC5Rsd9SMe79X+XjN0n6jcNfruslPdzMflfF/NH3R77vI1X0g78p/75YxcnyFZKur/y481YV+QKk4ketZ1XCsu+j4irfIu9x969K+qqZHVBxpf8pKq7qf6zsD/dV8aOPVIxt76xtY/MHoVZzaVPCFbWZcfd/lHR7ZQ7CC1T8Mhp6zY+Uv2I8w8yOVjGwvqRy8EnFgfuTm3+YWdOH8YdU/rphZt+k4gC8zsweruLA/R0Vc1u+dcFrr5P08MrfO8zs8eXjH5L04bK8y36lv0zSD5dzLh4n6R/d/eZy2TepuLwPTIpjNczM9prZyYFVPirpu81sWzkf4Wwdqr8PSfp5SR9UGQIl6atlnd+rZRv8haQfLB9X57n+maQf3fziaGYnm9mDlu3fAP5CRYKO+5RlaTv/60OSnmdmR5rZdhVXFj+6YL0PSnqOmd23vGr5fZJUXl37vJn9gHTvnLfQL9mMxXnr8sPBwh+fyitS7y3/3vzx6awFPz5tjh8nV5ZF/fgk6Vkqfhh4r5k9ecFqB21rwp/GH5/6vKJWXn36fRUhjv9OxZXE+3TZVlfufruKq05XqLj6+EfBF6zGJH1/pS13uPu1Desu6mOm4grh5usf6e7nl+v8a+2qvTSDH4Q4UZuncyT9ppl9SkUscps0yP9eRWjAr1YGoAer+DVpw4rJ8p9WcTAv8vuSjjCzq1WEDb2w/EXkLEnXlCExj1YRHlH3HhXzITZdJ+lFZnatpK9TMTE0xntV/EL0ORWD3n+pLHtS+T5ACtb6WDWzrzezGyX9rKRfNrMbzewB5RemR6iYmL9Q+ePLeSqyen1S0pXu/q5y8YdUXFX/YPnB/QWVJ48LxLbBi8t9vFrFnLfNcrxfRXjlX5bL3qFijtiovEjKcpmKuSj/U0XY2T8GX7TVpeVrP6liTssveBFGWn+fj6voL58s3+djlcXPl3SumX1SxVWBZwfej7F4ZvjxKazhx6fNk7IvlT+wPHfZdhb4iIofrU4o6/AHKsv+j7b+wHRY4hArkhEd4e7vlPTLkh4beK+v6ND4dp2KK6mPKP/ebO/rVFyh21k+/7zK6/9M0k9aeTnMzB6zYLsheyU9d/PHMDM73szqYeBV+f8g5AlMlOMf/9xdKuYzXF4+3inpmgHe44MqJzfzj3/86/Zv6GNVxQniK6bez9z+STq2/P9+KuaSPHbqMjWU80QVc+4mLwv/OrXflmNexdXr88vH1WQif6qtyUQ2k3pt04IEOpK+W9K/amtCkQeX619SbvPTkv6gXP98VRJ2qDmZyHkqfji4StL7VIQ51t/7v0r6z5X9+4yKkMdrVYTTRSUTUZHM6EYV89ruKB8/QMWFkRu0ODnTBZL+VsVV8TdU6vIidUsmskeHkol8g2rJRMrnX1hZ5zQV82s36/zpgff5fsUlE/k+HUom8geS3lw+f18V89euLttkM8nI8Sp+9KkmE6m27TWSdpaPn6dDiWGuVJmwSWWCp8prTlQxv3byY2aVf5uZZYAkmNlZKgdSFQfwo3vc9nYVE1D/tK9tAutqyGMV3ZjZW1QkEriPivCgX5+4SAuZ2bdL+jc/NMcJmFQZcvlGd39qeSWo7+8fj5b0o+7+s31tM2Vmdqy731leOXu1pM+6+ytHLsPPSPqyu79uzPftGydqAAAAWGv8+NSf8iTpHBWp9z8h6cfc/Z9HLsOPSHqTF1mBs8WJGgAAAJA4KzLyfk3t6Re4+9U9v8+rVaT7r3qVuy9M3Y/hcKIGAAAAAIkh6yMAAAAAJIYTNQAAAABIDCdqAAAAAJAYTtQAAAAAIDGcqAEAAABAYv5/89WlUqto7lIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x1080 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams['figure.figsize'] = (15, 15)\n",
    "fig = plt.figure()\n",
    "ax = [fig.add_subplot(131), fig.add_subplot(132), fig.add_subplot(133)]\n",
    "\n",
    "for i, (order, all_dofs_together, label) in enumerate([(0,False, \"non-zeros (p=0)\"),\n",
    "                                                    (1,False,\"non-zeros (p=1, low order + high order)\"),\n",
    "                                                    (1,True,\"non-zeros (p=1, all_dofs_together)\")]):\n",
    "    a = BilinearForm(L2(mesh,order=order,dgjumps=True,all_dofs_together=all_dofs_together))\n",
    "    a.Assemble()\n",
    "    ax[i].spy(sp.csr_matrix(a.mat.CSR()),markersize=3,precision=-1)\n",
    "    ax[i].set_xlabel(label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
